### "Allegiance, Ability, and Achievement in the American Civil War:
### Commander Traits and Battlefield Military Effectiveness"

### By Jeffrey B. Arnold, J. Tyson Chatagnier, and Gary E. Hollibaugh, Jr.

##### MISCELLANEOUS UTILITY FUNCTION CODE.
##### DO NOT CALL DIRECTLY.


# Bootstrapping function
mod.boot <- function(data, indices, model){
  d <- data[indices,]
  fit <- update(model, data = d)
  return(c(coefficients(fit)))
}

# Function to create bootstrapped covariance matrix
bootcovgen <- function(boot.out){
  boot.vcov <- matrix(NA, nrow = length(boot.out$t0),
                      ncol = length(boot.out$t0))
  
  for(i in 1:dim(boot.vcov)[1]){
    for(j in 1:dim(boot.vcov)[2]){
      boot.vcov[i,j] <- sum((boot.out$t[,i] - mean(boot.out$t[,i], na.rm = TRUE))*
                              (boot.out$t[,j] - mean(boot.out$t[,j], na.rm = TRUE)),
                            na.rm = TRUE)/(boot.out$R - 1)
    }
  }
  return(boot.vcov)
}


# Function to create texreg table for regression objects
extract.polr <- function(model, include.thresholds = FALSE, include.aic = TRUE,
                         include.bic = TRUE, include.loglik = TRUE,
                         include.deviance = FALSE,
                         include.nobs = TRUE, include.lrtest = TRUE, ...) {
  s <- summary(model, ...)
  lrtest.chi <- round(lmtest::lrtest(model)$C[2], 3)
  lrtest.pval <- round(lmtest::lrtest(model)$P[2], 3)

  tab <- s$coefficients
  zeta.names <- names(s$zeta)
  beta <- tab[!rownames(tab) %in% zeta.names, ]
  thresh <- tab[rownames(tab) %in% zeta.names, ]

  if (sum(!rownames(tab) %in% zeta.names) == 1) {
    beta <- t(beta)
    rownames(beta) <- rownames(tab)[!rownames(tab) %in% zeta.names]
  }
  if (sum(rownames(tab) %in% zeta.names) == 1) {
    thresh <- t(thresh)
    rownames(thresh) <- rownames(tab)[rownames(tab) %in% zeta.names]
  }

  threshold.names <- rownames(thresh)
  threshold.coef <- thresh[, 1]
  threshold.se <- thresh[, 2]
  threshold.zval <- thresh[, 1] / thresh[, 2]
  threshold.pval <- 2 * pnorm(abs(threshold.zval), lower.tail = FALSE)

  beta.names <- rownames(beta)
  beta.coef <- beta[, 1]
  beta.se <- beta[, 2]
  beta.zval <- beta[, 1] / beta[, 2]
  beta.pval <- 2 * pnorm(abs(beta.zval), lower.tail = FALSE)

  if (include.thresholds == TRUE) {
    names <- c(beta.names, threshold.names)
    coef <- c(beta.coef, threshold.coef)
    se <- c(beta.se, threshold.se)
    pval <- c(beta.pval, threshold.pval)
  } else {
    names <- beta.names
    coef <- beta.coef
    se <- beta.se
    pval <- beta.pval
  }

  n <- nobs(model)
  lik <- logLik(model)[1]
  aic <- AIC(model)
  bic <- BIC(model)
  dev <- deviance(model)
  lrt <- lrtest.chi
  lrp <- lrtest.pval

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.aic == TRUE) {
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.deviance == TRUE) {
    gof <- c(gof, dev)
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  if (include.lrtest == TRUE) {
    gof <- c(gof, lrt, lrp)
    gof.names <- c(gof.names, "Likelihood Ratio Test", "Likelihood Ratio p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }

  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }

  tr <- createTexreg(
    coef.names = names,
    coef = coef,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

setMethod("extract", signature = className("polr", "MASS"),
          definition = extract.polr)

# Function to get Quasi-AIC for quasibinomial models
x.quasibinomial <- function(...) {
  res <- quasibinomial(...)
  res$aic <- binomial(...)$aic
  res
}


# Function to get Quasi-AIC for quasiGLM objects
manual_QAIC <- function(x){QAIC(update(x, family = x.quasibinomial),
                                chat = deviance(update(x, family = x.quasibinomial))/
                                  df.residual(update(x, family = x.quasibinomial)))}


# Extension of the existing texreg functions for glm objects to use the quasi-AIC
extract.glm <- function(model, include.aic = FALSE, include.bic = FALSE,
                        include.loglik = FALSE, include.deviance = FALSE,
                        include.qaic = TRUE,
                        include.wald = TRUE,
                        include.nobs = TRUE, ...) {
  s <- summary(model, ...)

  coefficient.names <- rownames(s$coef)
  coefficients <- s$coef[, 1]
  standard.errors <- s$coef[, 2]
  significance <- s$coef[, 4]

  aic <- AIC(model)
  bic <- BIC(model)
  qaic <- manual_QAIC(model)
  waldstat <- lmtest::waldtest(model, test = "Chisq")$C[2]
  waldpval <- lmtest::waldtest(model, test = "Chisq")$P[2]
  lik <- logLik(model)[1]
  dev <- deviance(model)
  n <- nobs(model)

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.aic == TRUE) {
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.qaic == TRUE) {
    gof <- c(gof, qaic)
    gof.names <- c(gof.names, "Quasi-AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.deviance == TRUE) {
    gof <- c(gof, dev)
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.wald == TRUE) {
    gof <- c(gof, waldstat, waldpval)
    gof.names <- c(gof.names, "Wald Test", "Wald Test p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }

  tr <- createTexreg(
    coef.names = coefficient.names,
    coef = coefficients,
    se = standard.errors,
    pvalues = significance,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}


# Function to get the bootstrapped standard errors in a texreg object
extract.glm.boot <- function(model, boot, nullmod, include.aic = FALSE, include.bic = FALSE,
                             include.loglik = FALSE, include.deviance = TRUE,
                             include.qaic = TRUE,
                             include.wald = TRUE,
                             include.nobs = TRUE, ...) {
  s <- summary(model, ...)

  coefficient.names <- rownames(s$coef)
  coefficients <- s$coef[, 1]
  standard.errors <- sqrt(diag(bootcovgen(boot)))
  significance <- 2-2*pnorm(abs(coefficients/standard.errors))

  aic <- AIC(model)
  bic <- BIC(model)
  qaic <- manual_QAIC(model)
  waldstat <- lmtest::waldtest(model, nullmod, test = "Chisq")$C[2]
  waldpval <- lmtest::waldtest(model, nullmod, test = "Chisq")$P[2]
  lik <- logLik(model)[1]
  dev <- deviance(model)
  n <- nobs(model)

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.aic == TRUE) {
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.qaic == TRUE) {
    gof <- c(gof, qaic)
    gof.names <- c(gof.names, "Quasi-AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.deviance == TRUE) {
    gof <- c(gof, dev)
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.wald == TRUE) {
    gof <- c(gof, waldstat, waldpval)
    gof.names <- c(gof.names, "Wald Test", "Wald Test p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }

  tr <- createTexreg(
    coef.names = coefficient.names,
    coef = coefficients,
    se = standard.errors,
    pvalues = significance,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}


# Texreg extension for lrm objects (Design or rms package)
extract.lrm <- function(model, include.pseudors = FALSE, include.lr = TRUE,
                        include.aic = TRUE, include.bic = TRUE,
                        include.nobs = TRUE, include.pvalue = TRUE,
                        include.loglik = TRUE, ...) {
  attributes(model$coef)$names <- lapply(attributes(model$coef)$names,
                                         function(x) gsub(">=", " $\\\\geq$ ", x))
  coef.names <- attributes(model$coef)$names
  coef <- model$coef
  se <- sqrt(diag(model$var))
  p <- pnorm(abs(model$coef / sqrt(diag(model$var))),
             lower.tail = FALSE) * 2

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()

  if (include.pseudors == TRUE) {
    pseudors <- model$stats[10]  # extract pseudo R-squared
    gof <- c(gof, pseudors)
    gof.names <- c(gof.names, "Pseudo R$^2$")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  if (include.aic == TRUE) {
    aic <- AIC(model) # extract AIC
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  if (include.bic == TRUE) {
    bic <- BIC(model) # extract BIC
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  if (include.loglik == TRUE) {
    loglik <- as.numeric(logLik(model))  # extract LR
    gof <- c(gof, loglik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  if (include.lr == TRUE) {
    LR <- model$stats[3]  # extract LR
    gof <- c(gof, LR)
    gof.names <- c(gof.names, "Likelihood Ratio Test")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  if (include.pvalue == TRUE) {
    pvalue <- model$stats[5]  # extract LR
    gof <- c(gof, pvalue)
    gof.names <- c(gof.names, "Likelihood Ratio Test P-Value")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  if (include.nobs == TRUE) {
    n <- model$stats[1]  # extract number of observations
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }


  tr <- createTexreg(
    coef.names = coef.names,
    coef = coef,
    se = se,
    pvalues = p,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}




# Create function to read MPlus files
extract.mplus.model.boot <- function (model, model.boot, competence = TRUE,
                                      summaries = c("Observations", "CFI", "TLI", "RMSEA_Estimate"),
                                      clusters = FALSE, escape.latex = FALSE, ...) 
{ model.boot.params <- model.boot$parameters$stdyx.standardized
  if (competence == FALSE){model.boot.params <- subset(model.boot.params, 
                                                   !(paramHeader %in% c("COMPETEN.BY","COMPETEN.ON")))}
  else {model.boot.params <- subset(model.boot.params, !(paramHeader %in% c("CSALOYAL.BY","CSALOYAL.ON")))}
  if (summaries[1] != "none") {
    stopifnot(all(summaries %in% colnames(model$summaries)))
    knownsummaries <- c("Title", "Mplus.version", "AnalysisType", 
                        "DataType", "Filename", "Estimator", "Observations", 
                        "Parameters", "ChiSqM_Value", "ChiSqM_DF", "ChiSqM_PValue", 
                        "ChiSqM_ScalingCorrection", "ChiSqBaseline_Value", 
                        "ChiSqBaseline_DF", "ChiSqBaseline_PValue", "LL", 
                        "UnrestrictedLL", "LLCorrectionFactor", "UnrestrictedLLCorrectionFactor", 
                        "CFI", "TLI", "AIC", "BIC", "aBIC", "RMSEA_Estimate", 
                        "RMSEA_90CI_LB", "RMSEA_90CI_UB", "RMSEA_pLT05", 
                        "WRMR", "ObsRepChiSqDiff_95CI_LB", "ObsRepChiSqDiff_95CI_UB", 
                        "PostPred_PValue", "DIC", "pD", "SRMR.Within", "SRMR.Between", 
                        "Entropy", "AICC", "ChiSqM_Mean", "ChiSqM_SD", "ChiSqM_NumComputations", 
                        "LL_Mean", "LL_SD", "LL_NumComputations", "UnrestrictedLL_Mean", 
                        "UnrestrictedLL_SD", "UnrestrictedLL_NumComputations", 
                        "CFI_Mean", "CFI_SD", "CFI_NumComputations", "TLI_Mean", 
                        "TLI_SD", "TLI_NumComputations", "AIC_Mean", "AIC_SD", 
                        "AIC_NumComputations", "BIC_Mean", "BIC_SD", "BIC_NumComputations", 
                        "aBIC_Mean", "aBIC_SD", "aBIC_NumComputations", "RMSEA_Mean", 
                        "RMSEA_SD", "RMSEA_NumComputations", "SRMR_Mean", 
                        "SRMR_SD", "SRMR_NumComputations", "SRMR.Within_Mean", 
                        "SRMR.Within_SD", "SRMR.Within_NumComputations", 
                        "SRMR.Between_Mean", "SRMR.Between_SD", "SRMR.Between_NumComputations", 
                        "SRMR", "NumFactors")
    knownsummaries.decimals <- rep(TRUE, length(knownsummaries))
    names(knownsummaries) <- names(knownsummaries.decimals) <- knownsummaries
    knownsummaries.decimals[c("Title", "Mplus.version", "AnalysisType", 
                              "DataType", "Filename", "Observations", "Parameters", 
                              "ChiSqM_DF", "ChiSqBaseline_DF")] <- FALSE
    use.decimals <- rep(TRUE, length(summaries))
    use.decimals[which(summaries %in% knownsummaries)] <- knownsummaries.decimals[summaries[which(summaries %in% 
                                                                                                    knownsummaries)]]
    summary.values <- as.numeric(model$summaries[, summaries])
  }
  else {
    summaries <- character(0)
    summary.values <- numeric(0)
    use.decimals <- logical(0)
  }
  summaries <- c("Average Variance Explained (AVE)", 
                 "Squared Inter-trait Correlation", summaries[-1], 
                 "Number of Commanders", summaries[1])
  summaries[summaries == "CFI"] <- "Comparative Fit Index"
  summaries[summaries == "TLI"] <- "Tucker-Lewis Index"
  summaries[summaries == "RMSEA_Estimate"] <- "Root Mean Square Error of Approximation"
  summaries[summaries == "Observations"] <- "Number of Observations"
  if(competence == TRUE){ave <- subset(mplusboot$parameters$r2, param == "COMPETEN")$est}
  if(competence == FALSE){ave <- subset(mplusboot$parameters$r2, param == "CSALOYAL")$est}
  summary.values <- c(ave, model.boot.params$est[model.boot.params$paramHeader == "COMPETEN.WITH"]^2,
                      summary.values[-1], model$data_summary$overall$NClusters, summary.values[1])
  use.decimals <- c(TRUE, TRUE, use.decimals[-1], FALSE, use.decimals[1])
  model.boot.params$paramHeader[model.boot.params$paramHeader == "COMPETEN.ON"] <- "Formative"
  model.boot.params$paramHeader[model.boot.params$paramHeader == "COMPETEN.BY"] <- "Reflective"
  model.boot.params$paramHeader[model.boot.params$paramHeader == "CSALOYAL.ON"] <- "Formative"
  model.boot.params$paramHeader[model.boot.params$paramHeader == "CSALOYAL.BY"] <- "Reflective"
  model.boot.params$param[model.boot.params$param=="LOG_WINPCT"] <- "Log Cumulative Win Percentage"
  model.boot.params$param[model.boot.params$param=="COMMAND_EX"] <- "Command Experience"
  model.boot.params$param[model.boot.params$param=="TOPOVERALL"] <- "Top Grade Reached"
  model.boot.params$param[model.boot.params$param=="LOG_CUMOPP"] <- "Log Cumulative Opponent Casualty Percentage"
  model.boot.params$param[model.boot.params$param=="CONFED"] <- "Confederate Commander"
  model.boot.params$param[model.boot.params$param=="PREV_OFFIC"] <- "Confederate Office-holding Bias"
  model.boot.params$param[model.boot.params$param=="CONFED_PAR"] <- "Confederate Partisan Bias"
  model.boot.params$param[model.boot.params$param=="PEACECONVE"] <- "Peace Conference Attendance"
  model.boot.params$param[model.boot.params$param=="SOUTHBORN"] <- "Southern"
  model.boot.params$param[model.boot.params$param=="NUMBATTLES"] <- "Battle Experience"
  model.boot.params$param[model.boot.params$param=="USMA"] <- "United States Military Academy"
  model.boot.params$param[model.boot.params$param=="PASSMIDSHI"] <- "Passed Midshipman"
  model.boot.params$param[model.boot.params$param=="OTHERMA"] <- "Other Military Academy"
  model.boot.params$param[model.boot.params$param=="MEXWAR"] <- "Mexican War Experience"
  model.boot.params$param[model.boot.params$param=="TIME"] <- "Military Service"
  model.boot.params$param[model.boot.params$param=="PREWAROFFI"] <- "Prewar Elected Official"
  model.boot.params$param[model.boot.params$param=="BORDERBORN"] <- "Border State"
  model.boot.params$param[model.boot.params$param=="FOREIGNBOR"] <- "Foreign"
  model.boot.params$param[model.boot.params$param=="CSA_KINBIA"] <- "Confederate Bias in Weighted Coefficient of Relationship (Blood Relatives)"
  model.boot.params$param[model.boot.params$param=="CSA_KIN_IN"] <- "Confederate Bias in Weighted Coefficient of Relationship (In-Laws)"
  model.boot.params$param[model.boot.params$param=="PREWARPLAN"] <- "Planter"
  model.boot.params$param[model.boot.params$param=="NATIVE_LIN"] <- "Home-State Lincoln Vote"
  model.boot.params <- subset(model.boot.params, 
                   !(paramHeader %in% c("Thresholds", "Residual.Variances", 
                                        "Intercepts", "COMPETEN.WITH")))
  model.boot.params$headers <- paste0(model.boot.params$paramHeader,": ", model.boot.params$param)
  tr <- createTexreg(coef.names = as.character(model.boot.params$headers), 
                     coef = model.boot.params$est, se = model.boot.params$se, 
                     pvalues = model.boot.params$pval, 
                     gof.names = summaries, gof = summary.values, 
                     gof.decimal = use.decimals)
  return(tr)
}






